# import numpy as np
# import matplotlib.pyplot as plt
# import numpy as np
# from skimage import data,color
# img = data.coffee()
# img=color.rgb2gray(img)#直接读为灰度图像
# f = np.fft.fft2(img)
# fshift = np.fft.fftshift(f)

# #取绝对值：将复数变化成实数
# #取对数的目的为了将数据变化到0-255
# s1 = np.log(np.abs(fshift))
# """
# 巴特沃斯低通滤波器
# """
# def butterworthPassFilter(image, d, n):
#     f = np.fft.fft2(image)
#     fshift = np.fft.fftshift(f)

#     def make_transform_matrix(d):
#         transfor_matrix = np.zeros(image.shape)
#         center_point = tuple(map(lambda x: (x - 1) / 2, s1.shape))
#         for i in range(transfor_matrix.shape[0]):
#             for j in range(transfor_matrix.shape[1]):
#                 def cal_distance(pa, pb):
#                     from math import sqrt
#                     dis = sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
#                     return dis

#                 dis = cal_distance(center_point, (i, j))
#                 transfor_matrix[i, j] = 1 / (1 + (dis / d) ** (2*n))
#         return transfor_matrix
#     d_matrix = make_transform_matrix(d)
#     new_img = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift * d_matrix)))
#     return new_img

# plt.subplot(221)
# plt.axis("off")
# plt.title('Original')
# plt.imshow(img,cmap='gray')

# plt.subplot(222)
# plt.axis('off')
# plt.title('Butter 100 1')
# butter_100_1=butterworthPassFilter(img,100,1)
# plt.imshow(butter_100_1,cmap='gray')
# plt.subplot(223)
# plt.axis('off')
# plt.title('Butter 30 1')
# butter_30_1=butterworthPassFilter(img,30,1)
# plt.imshow(butter_30_1,cmap='gray')
# plt.subplot(224)
# plt.axis('off')
# plt.title('Butter 30 5')
# butter_30_5=butterworthPassFilter(img,30,5)
# plt.imshow(butter_30_5,cmap='gray')
# plt.show()


from skimage import data, color, util, io
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import filters

# 使用imread方法读入lena图片数据
img = io.imread('lena.bmp')  # 将咖啡图像换成lena图像
gray_img = color.rgb2gray(img)  # 将lena图像转换为灰度图像

# 为灰度图像添加高斯加性噪声
noisy_img = util.random_noise(gray_img, mode='gaussian')

# 在一个figure中显示图片
plt.figure(figsize=(20, 10))

# 显示原始灰度图像
plt.subplot(231), plt.imshow(gray_img, 'gray'), plt.title('原始灰度图像')

# 显示含噪图像
plt.subplot(232), plt.imshow(noisy_img, 'gray'), plt.title('含噪灰度图像')

# 巴特沃斯低通滤波器函数
def butterworthPassFilter(image, d, n):
    f = np.fft.fft2(image)
    fshift = np.fft.fftshift(f)

    def make_transform_matrix(d):
        transform_matrix = np.zeros(image.shape)
        center_point = tuple(map(lambda x: (x - 1) / 2, image.shape))
        for i in range(transform_matrix.shape[0]):
            for j in range(transform_matrix.shape[1]):
                def cal_distance(pa, pb):
                    return np.sqrt((pa[0] - pb[0])**2 + (pa[1] - pb[1])**2)
                dis = cal_distance(center_point, (i, j))
                transform_matrix[i, j] = 1 / (1 + (dis / d)**(2*n))
        return transform_matrix
    d_matrix = make_transform_matrix(d)
    new_img = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift * d_matrix)))
    return new_img

# 设置图像尺寸和截止频率
img_shape = noisy_img.shape
# 中文显示工具函数
def set_ch():
    from pylab import mpl
    mpl.rcParams['font.sans-serif'] = ['FangSong']
    mpl.rcParams['axes.unicode_minus'] = False
set_ch()

# 应用巴特沃斯滤波器，并显示结果
# 不同截止频率和阶数的巴特沃斯滤波器
butter_100_1 = butterworthPassFilter(noisy_img, 100, 1)
plt.subplot(233), plt.imshow(butter_100_1, 'gray'), plt.title('巴特 100 1')

butter_30_1 = butterworthPassFilter(noisy_img, 30, 1)
plt.subplot(234), plt.imshow(butter_30_1, 'gray'), plt.title('巴特 30 1')

butter_30_5 = butterworthPassFilter(noisy_img, 30, 5)
plt.subplot(235), plt.imshow(butter_30_5, 'gray'), plt.title('巴特 30 5')

# 使用 scipy.ndimage.filters.gaussian_filter 进行比较
gaussian_denoised = filters.gaussian_filter(noisy_img, sigma=1)
plt.subplot(236), plt.imshow(gaussian_denoised, 'gray'), plt.title('高斯滤波去噪')

plt.tight_layout()
plt.show()
